Surface Solitons in Three Dimensions 
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We study localized modes on the surface of a three-dimensional dynamical lattice. The stability 
of these structures on the surface is investigated and compared to that in the bulk of the lattice. 
Typically, the surface makes the stability region larger, an extreme example of that being the three- 
site "horseshoe" -shaped structure, which is always unstable in the bulk, while at the surface it is 
stable near the anti-continuum limit. We also examine effects of the surface on lattice vortices. For 
the vortex placed parallel to the surface this increased stability region feature is also observed, while 
the vortex cannot exist in a state normal to the surface. More sophisticated localized dynamical 
structures, such as five-site horseshoes and pyramids, are also considered. 



INTRODUCTION 



Surface waves have been a subject of interest in a vari- 
ety of contexts, including surface plasmons in conductors 
[if and optical solitons in waveguide arrays Q in physics, 
surface waves in isotropic magnetic gels |3[ in chemistry, 
water waves in the ocean in geophysical hydrodynamics, 
and so on. Quite often, features exhibited by such wave 
modes have no analog in the corresponding bulk media, 
which makes their study especially relevant. In particu- 
lar, a great deal of interest has been drawn to nonlinear 
surface waves in optics. It was shown theoretically [i| 
and observed experimentally [f| that discrete localized 
nonlinear waves can be supported at the edge of a semi- 
infinite array of nonlinear optical waveguide arrays. Such 
solitary waves were predicted to exist not only in self- 
focusing media (as in the above-mentioned works), but 
also between uniform and sclf-dcfocusing media 0, @], 
or between self-focusing and self-defocusing media (e.g. 
in Q). They have been subsequently observed in me- 
dia with quadratic Q and photorefractive [^, nonlin- 
earities. In the two-dimensional (2D) geometry, stable 
topological solitons have been predicted in a saturable 
medium [Til ], which constitute generalizations to lattice 
vortex solitons predicted in Ref. [l2| • Quasi-discrete vor- 
tex solitons have been experimentally observed in a self- 
focusing bulk photorefractive medium [l3| . Theoretical 
predictions for a variety of species of discrete 2D sur- 

l&k , corner modes Wl UJt , 



or a variety ot 
face solitons M, (Tfl Gil, 
as well as surface breathers 



17| , were reported. Subse- 



quent work has resulted in the experimental observation 
of 2D surface solitons. both fundamental ones and multi- 



pulse states, in photorefractive media [l9(, as well as in 
asymmetric waveguide arrays written in fused silica [2(| • 
Recently, surface solitons in more complex settings, such 
as chirped optical lattices in ID and 2D [2l|, [22j , at inter- 
faces between photonic crystals and metamaterials [23j |. 
and in the case of nonlocal nonlinearity [24|, Hfl, have 
emerged. 

Nearly all these efforts have been aimed at the study 
of surface solitons in ID and 2D geometries. The only 
3D setting examined thus far assumed a truncated bun- 
dle of fiber-like waveguides, incorporating the temporal 
dynamics in longitudinal direction to produce 3D "sur- 
face light bullets" in Ref. [26| (the respective 2D surface 
structures were examined in Ref. [27j). 

Our aim in the present work is to extend the analysis to 
surface solitons in genuine 3D lattices. Our setup is rel- 
evant to a variety of applications including, e.g. , crystals 
built of microresonators trapping photons [281 ] . polari- 
tons [2^], or Bose-Einstein condensates in the vicinity of 
an edge of a strong 3D optical lattice [10, HH . In partic- 
ular, we report results for discrete solitons at the surface 
of a 3D lattice, i.e., 3D localized states that are similar 
to relevant objects studied in the 2D setting of Ref. [HI]. 
Thus, we will study localized states such as dipoles and 
"horseshoes" abutting on a set of three lattice sites, but 
also states that are specific to the 3D lattice. A vari- 
ety of species of such solitons is examined below, and 
their stability on the surface is compared to that in the 
bulk. Some localized structures, such as dipoles, may 
be placed either normal or parallel to the surface. We 
demonstrate that, typically, the enhanced contact with 
the surface increases the stability region of the struc- 
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ture. Physically, this conclusion may be understood by 
the fact that the surface reduces the local interactions to 
fewer neighbors, rendering the system "more discrete", 
hence more stable (by pushing the medium further away 
from the continuum limit, where all solitons would be 
unstable against the collapse). This effect is remarkable, 
e.g., for the three-site horseshoes which are never sta- 
ble in the bulk, but get stabilized in the presence of the 
surface. However, the surface may also have an adverse 
effect, inhibiting the existence of a particular mode. The 
latter trend is exemplified by discrete vortices, which, if 
placed parallel to the surface, feature enhanced stability 
as compared to the bulk-mode counterpart, but cannot 
exist with the orientation perpendicular to the surface. 
Surface-induced effects of a different kind, which are less 
specific to discrete systems, are induced by the inter- 
action of a particular localized mode with its fictitious 
"mirror image" . In terms of lattice models, the approach 
based on the analysis of the interaction of a real mode 
with its image was proposed in Ref. [32j |. 

To formulate the model, we introduce unit vectors 
ei = (1,0,0), e 2 = (0,1,0), and e 3 = (0,0,1) and de- 
fine lattice sites by n = n j e j with integer nj. We 
assume that the lattice occupies a semi-infinite space, 
n 3 > 1; an d its dynamics obeys the discrete nonlinear 
Schrodingcr (DNLS) equation in its usual form, 

i4> n +e^4> n + a\4> n \ 2 (j) n = {). (1) 

Here </> n is a complex discrete field, e is the coupling con- 
stant, n stands for the time derivative, the parameter 
a = ±1 determines the sign of the nonlincarity (focusing 
or defocusing respectively), and A0 n is the 3D discrete 
Laplacian: 

3 

n+ej + 0n-e 3 — 2s0 n ) , (2) 

i=i 

for 713 > 2, while for 713 = 1 the term with subscript index 
n — e3 is to be dropped (note that e3 is the direction 
normal to the surface). 

It is interesting to point out here that an approach to- 
wards understanding the dynamics of Eq. ([1]) in the vicin- 
ity of the surface can be based on the above-mentioned 
concept of the fictitious mirror image, formally extends 
the range of 713 up to 713 = — 00, supplementing the equa- 
tion with the anti-symmetry condition, 

4 > n 1 ,n2. 1 —n 3 = 4 > ni,n2n 3 - (3) 

Indeed, this condition implies 4> ni ,n 2 ,o = 0, which is 
equivalent to the summation restriction in Eq. ([2]) as de- 
fined above. 

To confine the analysis to localized solitary wave 
modes, we impose zero boundary conditions, n — > 
at 71^2 — > ±oo and 713 — * 00. Additionally, s = ±1 
in Eq. ^ — this parameter is introduced for convenience 
(see Sec. lIII Bj) and can be freely rescaled using the trans- 
formation <f> —* <f> e lvt for an appropriate choice of v and 



time rescaling. Stationary solutions to Eq. |T]) will be 
sought for as </>„ = exp (iAt)u n , where A is the frequency 
and the lattice field u n obeys the equation 

(A-a\u n \ 2 )u a ~-eAu a = 0. (4) 

Our presentation is structured as follows. The follow- 
ing section recapitulates the necessary background for the 
prediction of the existence and stability of lattice solitons. 
In section III, we report a bifurcation analysis for various 
surface states, treated as functions of coupling constant 
e, with emphasis on the comparison with bulk counter- 
parts of these states. Section IV reports the study of the 
evolution of unstable surface states. Finally, section V 
summarizes our findings and presents our conclusions. 

II. THE THEORETICAL BACKGROUND 

First, we outline some general properties of the model. 
Equation ([1]) conserves two dynamical invariants, namely 
the norm N , 

00 

^V= E l^l 2 ' ( 5 ) 

"3 = 1 

"•1,2=— 00 

and the Hamiltonian H , 

00/3 \ 

H= Y, UE^n(^n + e J -^„)+C.C.]+||0 n | 4 
"3 = 1 \ 3=1 I 

m 2=— 00 

(6) 

where the asterisk stands for complex conjugation. Sta- 
tionary solutions to Eq. (01 with a = ±1 are connected 
by the staggering transformation [l7l . [33| : if u n is a solu- 
tion for some A and a = +1, then (— l)" 1+,l2+ " 3 w n is a 
solution for A = 12s — A and a = — 1. Consequently, it is 
sufficient to perform the analysis of stationary solutions, 
including their stability, for a single sign of the nonlin- 
carity; thus, below we will fix a = +1 (corresponding to 
the case of onsite self-attraction). 

Solutions to Eq. (H} in half-space 713 > 1, subject to 
boundary condition n = for 713 = 0, as defined above, 
may be continued anti-symmetrically for the entire 3D 
space by setting U n = u n for 773 > 1 and U n = —u n for 
?7 3 < —1. Then, according to results of Ref. [34|], this 
leads to an immediate conclusion, namely that there ex- 
ists a minimum norm N m in necessary for the existence 
of localized surface states in the present model. In other 
words, no surface modes survive in the limit of N —> 0. 
In this connection, it is relevant to note that numerical 
findings that will be presented below were obtained, of 
course, for finite cubic lattices where, strictly speaking, 
there is no lower limit for N necessary for the existence 
of localized modes [l?}- At this point, we have to specify 
that speaking about localized modes in a finite lattice 
we understand solutions which are localized on a num- 
ber of cites much smaller than the total number of sites 
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in the chosen direction used for numerical simulations. 
Next we recall that generally speaking, there exist sev- 
eral branches of the nonlinear localized modes, i.e. for 
a given e one can find localized excitations at different 
values of the norm N. Using the natural terminology 
we refer to higher/lower branches speaking about solu- 
tions with larger/smaller norm. In this classification the 
surface modes we are dealing with correspond to higher 
branches of the solutions of the respective finite lattices, 
i.e., their norm cannot be made arbitrarily small (sec also 
the relevant discussion below in Section UlI Bp . 

To find solution families, we start with the anti- 
continuum (AC) limit, e = [35|. In this limit, the 
lattice field is assumed to take nonzero values only at 
a few ("excited") sites, which determines the profile of 
the configuration to be seeded. The continuation of the 
structure to e > is determined by the Lyapunov's re- 
duction theorem (36[. More specifically, the solution is 
expanded as a power series in e, the solvability condition 
at each order being that the respective projection to the 
kernel generated by the previous order does not give rise 
to secular terms |35f . 

The linear stability is then studied, starting from the 
usual form of the perturbed solution, 



iAI 



(u n + Sa n e 



XI 



(7) 



where S is a formal small parameter, and A is a stabil- 
ity eigenvalue associated with eigenvector ip — {a n ,6*}. 
Substituting this into Eq. (fT]) yields the linearized system, 



iXa-a = — eAc 



2|u r 



ifb* 



-iXbl = -eAb* n + Ab* n -2\u n \ 2 b* n ~u* n 2 a a , 

which can be cast in the form 
/ W (i,a) 

I H (2,l) H (2,2) 




(8) 



where A and B are vectors composed by elements a n and 
&*, respectively, while the matrices Ti^^ (p,q £ {1,2}) 
are given by, 



= = <W (A + 6 S£ - 2\ Ul 



rl n,n' ri n,n> ~ "n,n"* n t . 



(9) 



An underlying stationary solution is (spectrally) unstable 
if there exists a solution to Eq. © with Rc(A) > 0. Oth- 
erwise, the stationary solution is classified as a spectrally 
stable one. As explained in Ref. [I?) , the Jacobian of the 
above mentioned solvability conditions is intimately con- 
nected to the full eigenvalue problem. More specifically, 
if the eigenvalues 7 of the M x M eigenvalue problem of 
the Jacobian (where M is the number of excited sites at 



the AC limit), then the near-zero eigenvalues of the full 
stability problem can be predicted to be A = ^f2rfe p ' 2 , 
where p is the number of lattice sites that separate the 
adjacent excited nodes of the configuration at the AC 
limit. 



III. THE BIFURCATION ANALYSIS 
A. Existence and stability of surface structures 

In this section we study, by means of numerical meth- 
ods, the existence and stability of various 3D configura- 
tions and compare the results to the corresponding ana- 
lytical predictions. These configurations are obtained by 
starting from the AC limit (e = 0), and are continued to 
£ > 0, using fixed-point iterations. For all the numerical 
results presented in this work, we fix the normalization 
A = 1 [see Eq. (j3}], and use a lattice of size 13 x 13 x 13, 
unless stated otherwise. Also, for the presentation of the 
numerical results, we replace the triplet (711,712,713) by 
(l,n,m), i.e., the surface corresponds to m = 1. 

We start by examining dipoles aligned parallel or nor- 
mal to the surface. Panel (a) in Fig. [T] shows the norm 
of such states versus coupling constant e, while panel (b) 
depicts the imaginary part of the stability eigenvalues 
for the bulk dipolc, produced by the theory outlined in 
the previous section [dashed (black) lines], and by the 
numerical computations [solid (blue) lines]. It is worth 
mentioning that, for all the different configurations that 
we report in this manuscript, we display the imaginary 
part of the stability eigenvalue only for the bulk mode 
since the difference between the curves for the different 
variants (bulk, parallel or normal to the surface) is min- 
imal. It should be noted however that the contact with 
the surface may produce higher order (smaller) eigenval- 
ues that are not present in its bulk counterpart (results 
not shown here). The theoretical prediction for the sta- 
bility eigenvalues is A = ±2y / S, which, as expected, is 
the same as in an out-of-phase (twisted) ID mode an- 
alyzed in Ref. [3^ ]. since the structure is nearly one- 
dimensional, along the line connecting the two excited 
sites. Panel (c) in Fig. [T] compares the largest instability 
growth rate as a function of e for the bulk dipoles [dash- 
dotted (green) line] and those oriented normally and par- 
allel to the surface [dashed (red) and solid (blue) lines, 
respectively] . It is seen that the stability interval of the 
dipoles increases as its contact with the surface strength- 
ens, in accordance with the arguments presented above. 
In the case of the bulk dipole, the instability occurs for 
values of the coupling constant in between £0 = 0.114 
and £1 = 0.115. From now on, when reporting computed 
instability thresholds, we will use the lower bound for e 
(e.g., £0 in the above example) with the understanding 
that we always used the same resolution in e. For the 
dipolc set normally to the surface, we observe the onset 
of instability at e = 0.117, while for the parallel-oriented 
one at e = 0.120. In panels (e)-(h) of Fig.Q]wc also depict 




FIG. 1: (Color Online) Results for the dipoles oriented paral- 
lel and normal to the surface, (a) Norm N versus the lattice 
coupling constant, e. (b) Imaginary part of the linear stability 
eigenvalue: solid (blue) and dashed (black) lines correspond, 
respectively, to numerically found and analytically predicted 
forms, (c) Real part of the critical (in)stability eigenvalue: 
the dashed (red) and solid (blue) lines depict the normal- and 
parallel-oriented dipoles, respectively, while the dash-dotted 
(green) line corresponds to the bulk dipole. (d) (In)stability 
eigenvalue for the parallel surface dipole placed at distances 
from the surface starting from zero and up to five lattice 
periods away (curves right to left), (e)-(g) Configurations 
and (f)-(h) respective spectral stability planes just above the 
instability threshold. The level contours, corresponding to 
Re(«;,n,m) = i0.5 max {?i;. n , m } are shown, respectively, in 
dark gray (blue) and gray (red). The instability thresholds 
for the dipoles oriented parallel and normally to the surface 
are, respectively, e = 0.117 and e = 0.120. For comparison, 
the threshold for the bulk dipole is e = 0.114. 



the shapes of the normal and parallel dipoles, just below 
the instability threshold, along with their corresponding 
spectral stability planes. 

The stabilizing effects exerted by the surface depend, 
in a great measure, on the distance of the configuration 
from the surface, namely, the further away the config- 
uration from the surface, the lesser the effect is. This 
property is clearly seen in panel (d) of Fig. [TJ where we 
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FIG. 2: (Color Online) The stability of the three-site "horse- 
shoe". Panels are similar to those in Fig. [T] Panel (c) 
compares the critical stability eigenvalue, as a function of 
the lattice coupling, e, for the surface and bulk horseshoes 
[solid (blue) and dashed-dotted (green) lines, respectively]. 
The bulk horseshoe is always unstable (due to a purely real, 
higher-order eigenvalue) , while the corresponding surface con- 
figurations have a stability region (the corresponding eigen- 
value becomes imaginary in this case). Panels (d)-(e) corre- 
spond to the surface horseshoe just above the stability thresh- 
old of e = 0.239. 



plot the (in)stability eigenvalue as a function of the cou- 
pling for several values of the separation of the parallel 
dipole from the surface. The curves, from right to left, 
depict the results for the dipole set at the distance of 
0, 1, 5 sites away from the surface (0 sites refers to 
the dipole sitting on the surface). As the panel demon- 
strates, the stability interval is reduced as the dipole is 
pulled away from the surface, converging towards a bulk 
dipole. 

Let us now consider the "horseshoe" configurations, 
for which the presence of the surface is critical to their 
stability. In Fig. [2] we depict the properties of a three- 
site horseshoe, which actually is a truncated version of a 
quadrupole, cf. the 2D situation [TJ]. As before, panel 
(a) in Fig. [2] shows the norm versus e, while panels (b) 
and (c) compare the stability of the bulk horseshoe (the 
dash-dotted line) and ones built near the surface (the 
solid line). While the bulk horseshoes are always unsta- 
ble, similar to their 2D counterparts [3], the ones placed 
near the surface are stable at small e, destabilizing at 
e = 0.239. Panels (d)-(e) in Fig. [2] show the configura- 
tion for the coupling just below the instability threshold, 
along with the respective spectral plane. The analytical 
expressions for stable eigenvalues are 
A = 0(e 2 ), cf. the expressions obtained in Ref. [3] for 
the 2D horseshoes. 
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FIG. 3: (Color Online) The stability for the five-site horseshoe 
at the surface. Panels are identical to those in Fig. [2] In this 
case, the stability threshold is at e — 0.211, while for the 
bulk 5-site horseshoe it is e — 0.205. Panels (d)-(e) depict 
the configuration and the respective linear stability spectrum 
just above the critical point of s = 0.211. 



Figure [3] illustrates the same features as before but for 
the five-site horseshoe. Unlike its three-site cousin, the 
bulk five-site horseshoe is stable up to a critical value 
of the coupling, e = 0.205, while the surface variant 
has it stability region e < 0.211. The eigenvalues of 
the linearization in this case can be computed similar 
to those for the three-site modes |14j . as outlined above 
(cf. also Ref. [HI), which eventually yields A = 3.8042ei, 
A = 2.8284ei, A = 2.3511ei, A = C(e 2 ), and A = 0, 
in good agreement with the corresponding numerical re- 
sults, as shown in panel (b) Fig. [3] 

Next we consider the quadrupole configuration, see 
Fig- SI The surface again exerts a stabilizing effect, albeit 
a weaker one, when the quadrupole is placed normally 
and parallel to the surface. In the bulk, the quadrupole 
loses stability at e = 0.068, while the normal and paral- 
lel surface quadrupolcs have stability thresholds, respec- 
tively, at e = 0.070 and e = 0.071. The analytical ap- 
proximation for the stability eigenvalues in this case are 
A = y/Bei (a double eigenvalue), A = 2^/ei, and a zero 
eigenvalue, which accurately capture the numerical find- 
ings depicted in panel (b) of Fig. [4] 

In Fig. [5] we present the results for four-site vor- 
tices. This configuration, in contrast to the previous 
ones, is described by a complex solution. In the AC 
limit, the vortex occupies the same excited sites as 
the above-mentioned quadrupole, but the phase profile, 
|0, 7r/2, 7t, 37r/2}, emulates that of the vortex of charge 1 
[121 , [351 ] . The bulk four-site vortex (which was discussed 
in Ref. (38|) loses its stability at e = 0.438, while the 
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FIG. 4: (Color Online) The stability of quadrupole modes. 
The layout is similar to that in Fig. [3] In panel (c), due to 
the close proximity of the thresholds, the close-up is shown 
for the critical stability eigenvalue versus the lattice coupling 
constant, e, for the parallel and normal surface modes, and 
the bulk one [solid (blue) and dashed (red) lines, and the 
dash-dotted (green) line, respectively]. The threshold for the 
bulk mode is e = 0.068, while for the normal and parallel 
quadrupoles it is, respectively, s = 0.070 and e = 0.071. As 
before, panels (d) and (e) show the configuration just above 
the instability threshold along with its corresponding spectral- 
stability plane. 



vortex parallel to the surface features an extended sta- 
bility region, up to e = 0.505. However, the surface in 
this case prohibits the existence of a vortex that would be 
oriented normally to the surface layer, similarly to what 
was found for 2D lattice vortices [14 ] . 

The simplest explanation for the complete absence of 
the solution normal to the surface, compared with that 
of an existing vortex waveform parallel to the surface can 
arguably be traced in the interaction of such vortices in 
the half-space with their fictitious image (if the domain 
is equivalently extended to the full space). In the case of 
the vortex parallel to the surface, the situation is tanta- 
mount to the vortex cube structures examined in (39l.[40|. 
for which it was established in [|(| that the persistence 
conditions are satisfied (and, in fact, that such structures 
consisting of two out-of-phase vortices should be linearly 
stable close to the AC limit). On the other hand, for 
the case normal to the surface, by examining the rele- 
vant interactions it can be observed (at an ap pro priatel y 
high order) that the persistence conditions of [35|, [33, |4(| 
can not be satisfied and hence the structure can not be 
continued beyond the AC limit. That is why the struc- 
ture can never be observed to exist irrespectively of the 
smallncss of e. 

The next species of stationary lattice solutions is a 
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FIG. 5: (Color Online) The stability of the four-site vortex in 
the grid of size llx llx 11. The dash-dotted and solid lines 
show the bulk vortex and the one parallel to the surface, re- 
spectively. The layout is similar to that of the above figures. 
Instability in the bulk occurs at e = 0.438, and in the parallel 
surface vortex at e = 0.505. The vortex cannot exist with the 
orientation normal to the surface. Panels (d) and (e) show 
the parallel surface vortex just above the instability thresh- 
old of e = 0.485. As in the previous figures, the level con- 
tours corresponding to Ke(ui t n ym ) = ±0.5 max {«;,„,,„} are 
shown, respectively, in dark gray (blue) and gray (red), while 
the complementary level contours, defined as Im(Wi ;Tl , m ) = 
±0.5 max {ti; jrl . m }, are shown by light gray (green) and very 
light gray (yellow) hues, respectively. 



pyramid-shaped structure, with characteristics displayed 
in Fig. [51 whose base is a rhombus composed of four 
sites. The remaining out-of-plane vertex site must have 
phase or ir, since the phase values ir/2 and 3ir/2 at this 
site do not produce a solution. The full set of pyramids 
(bulk, normal, parallel — see panels (d)-(f) of Fig. [6]) is 
completely unstable, as seen in panel (c) of Fig. [fJl the 
surface producing no stabilizing effect on it. This strong 
instability actually arises at the lowest order in the an- 
alytical eigenvalue calculations, which yield A = 2\Z5ei, 
A = 2\/2ei, A = 2e, A = 0, and A = 0(e 2 ), once again 
in very good agreement with the full numerical results of 
Fig. E 



B. Small-amplitude modes in a finite lattice 

Since our numerical investigation of the surface modes 
uses a finite lattice, which allow the existence of small- 
amplitude modes (ones with the zero threshold in terms 
of the norm — cf. discussion in Sec. [TTJ) , here we briefly 
consider the modes in a finite lattice having the small- 
amplitude limit. Our aim is to show that these modes 
belong to lower branches, as compared with the "nor- 
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FIG. 6: (Color Online) The instability of pyramid-shaped 
structures. This configuration abuts on the base in the form of 
a rhombus, and includes the out-of-plane site with zero phase. 
Three variants of this configuration are displayed in panels 
(d)-(f): bulk, normal and parallel to the surface, respectively. 
The stability of the three different variants of the pyramid is 
essentially identical, all three of them being unstable. 



mal" surface modes considered above. To this end, we 
concentrate on the lattice of size Mx Mx M lattice, sub- 
ject to the zero boundary conditions, which imply that 
discrete Laplacian ^ is modified at surfaces rij = 1 and 
rij = M (j = 1, 3) by setting the fields at sites n — 
and n + ej , respectively, equal to zero. For the sake of 
definitcness, we fix here s = — 1 in Eq. @. 

To determine the norm TV of small-amplitude modes 
we follow Ref. (l7j . and look for a solution to Eq. ((4]) 
with the amplitude it n and coupling constant e being 
represented as series 



U n = euo.n ± e 2 «2,n ± C(e 3 ), 



s + e 2 e 2 + O(e i ), 



(10) 



in powers of small parameter e = y/%N/ (M + l) 3 <C 1, 
which vanishes in the limit of the infinite lattice (M — > 
oo); in other words, small e characterizes the "largeness" 
of the lattice. We focus on real solutions here. 

Substituting expansions (jXOj) into Eq. (j4]) and gather- 
ing terms of the same order in e, we rewrite Eq. ((¥]) in 
the form of a set of equations: 



Art,- „ - e Att,- n = F, n . 



(11) 



Here F 0i „ = 0, F 2i „ = A(e 2 /e )u . n + (u ,„) , hence 
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Eq. (TTTj) with ,j = gives rise to a linear eigenmodc, 



3=1 



M+ 1 / ' 



(12) 



with the respective approximation for the lattice coupling 
constant, 



'0 



A 



6 + 2£< 

3=1 



71/ + 1 



(13) 



parameterized by vector m = (m l5 m 2 , m 3 ). At the same 
time, considering the solvability conditions for j = 2, 
which amounts to demanding the orthogonality of i<2,n 
and i*o. n , we obtain corrections to the coupling constants, 



c 2 = .( m ) 3 



64A 



n(3+* 



m j ,(M+l)/2) 



(14) 



It follows from Eq. (j 14[) that each of the linear modes 
(|12|) is uniquely extended into a small-amplitude nonlin- 
ear one. These modes are characterized by the linear 
dependence of the norm on coupling constant e: 



8A(M + l) 3 



ei m) - e 



"O^ Ilj=l ( 3 + $m j ,(M+l)/2) 



(15) 



From Eq. ([15)) it follows that each mode, parameter- 
ized by vector m, exists when e belongs to the inter- 
val < e < e Q m ' > . The validity of approximation (fl~5|) is 
corroborated by the coincidence of analytical and numer- 
ical results in the vicinity of £q (as shown in Fig. [7]), 
where these modes reaches their small- amplitude limit. 
Such a property of these modes differs considerably from 
the case of the surface modes which do not possess the 
small-amplitude limit and require some minimal value of 
the norm (for the normal dipole, depicted in Fig. [7] by 
dash-dotted line, the minimal norm is r* 1.262). Panel 
(b) in Fig. [7] shows that only the mode, parameterized 
by vector m = (1, 1, 1), is stable for e close to £ m \ while 
other modes are completely unstable. 



IV. DYNAMICS 

In this section we examine the nonlinear evolution of 
the various configurations, displaying the results in a set 
of figures (see Figs. [8l fT2|) . In each case, the evolution is 
initiated at a value of the coupling e taken beyond the 
instability threshold, and an initial small random per- 
turbation is applied in order to expedite the onset of the 
instability. 

All the figures display the evolution of the instability 
at six different moments of time, starting at t — 0, and 
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FIG. 7: (Color Online) Low-amplitude modes in a finite grid 
of size 9x 9x 9 with A = 1.0. (a) Norm N versus coupling con- 
stant e for several modes whose low-amplitude limit is param- 
eterized by vectors m, as calculated numerically and predicted 
by approximation (|15|l (solid and dashed lines, respectively). 
For comparison, dash-dotted line depicts the norm for sur- 
face normal dipole. (b) Real part of the critical (in)stability 
eigenvalue, calculated numerically. 



ending at a time well beyond the point at which the insta- 
bility manifests itself. All configurations that were pre- 
dicted above to be unstable through nonzero real parts 
of the (in)stability eigenvalue A indeed exhibit instabil- 
ity dynamics, which eventually results in a transition to 
a different configuration. In the case of the dipoles and 
horseshoes, Figs. [8rH0l show a spontaneous transition to 
monopolc patterns, i.e., ones centered around a single ex- 
cited site. On the other hand, in the case of the vortices 
and pyramids shown in Figs. 11111121 a few sites may re- 
main essentially excited at the end of the evolution. The 
monopole is, obviously, the most robust dynamical state 
in the lattice system, with the widest stability interval, in 
comparison with other discrete structures. This simplest 
state becomes unstable, for given A, only at values of the 
coupling constant e m A |38| . Another structure with a 
relatively wide stability region is the dipole (the more sta- 
ble the wider the distance between its constituent sites 
[1^]), consonant with the observation that some of the 
structures (especially ones with a large number of ex- 
cited sites, such as vortices and pyramids) dynamically 
transform into dipoles. 

Generally speaking, the exact scenario of the nonlin- 
ear evolution and the finally established state depend on 
details of the initial perturbation. In the figures, each 
configuration is shown by iso-level contours of distinct 
hues. In particular, dark gray (blue) and gray (red) are 
iso-contours of the real part of the solutions, while the 
light gray (green) and very light gray (yellow) colors de- 
pict the imaginary part of the same solutions. 

A case that needs further consideration is that of the 
three-site horseshoe. As observed from the stability anal- 
ysis presented in Fig. [21 this horseshoe in the bulk gives 
rise to a small unstable purely real eigenvalue for all val- 
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FIG. 8: (Color Online) The evolution of unstable dipoles: (a) 
a bulk dipole; (b) and (c) dipoles placed parallel and normally 
to the surface, respectively. In all the cases, the dipole is 
subject to oscillatory instability, which is responsible for the 
eventual concentration of most of the norm at a single site 
(i.e., the transition to a monopole). Parameters are A = 1, 
e — 0.2, the lattice has a size of 13 x 13 x 13, and times are 
indicated in the panels. All iso-contour plots are defined as 
Re(«i,n,m) = ±0.75 = Im(u; in , m ), and the initial configura- 
tions were perturbed with random noise of amplitude 0.01. 
The coding for the iso-contours is as follows: dark gray (blue) 
and gray (red) colors pertain to iso-contours of the real part of 
the solutions, while the light gray (green) and very light gray 
(yellow) colors correspond to the iso-contours of the imagi- 
nary part. 



ues of £, see the lower green dashed-dotted curve in panel 
(c) of the figure. Despite the presence of this eigenvalue, 
the evolution of the unstable bulk three-site horseshoe is 
predominantly driven by the unstable complex eigenval- 
ues, if any (in fact, for e > 0.226, see the dashed-dotted 
(green) line of Fig. [3J(c)). A careful analysis of the insta- 
bility corresponding to the small purely real eigenvalue 
for £ < 0.226 (i.e., before the complex eigenvalues be- 
come unstable) reveals that the corresponding dynamics 
amounts to an extremely weak exchange of the norm be- 
tween the two in-phase excited sites (see Fig. [5]). The 




FIG. 9: (Color Online) The evolution of the unstable three- 
site horseshoes: (a) bulk three-site horseshoe and (b) the 
horseshoe oriented normally to the surface. In both cases, 
the unstable horseshoe is subject to an oscillatory instability, 
which leads to the eventual concentration of most of the norm 
in a single-site structure. The iso-contours and parameters are 
the same as in Fig. [H] except that e = 0.3. 




FIG. 10: (Color Online) The evolution of unstable five-site 
horseshoes: (a) the bulk horseshoe, and (b) the five-site horse- 
shoe oriented normally to the surface. In both cases, the un- 
stable horseshoe is subject to an oscillatory instability, which 
triggers the transition to a monopole. The iso-contours and 
parameters are the same as in Fig. [9] 
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FIG. 11: (Color Online) The evolution of unstable vortices: 
(a) the bulk vortex for e = 0.3 and (b) the vortex parallel to 
the surface, for e — 0.6 and A = 1. The iso-contour plots are 
defined by Re(ui jn , m ) = ±1 = Im(Mj, n , m ). 



norm exchange is driven by the corresponding unsta- 
ble eigcnfunction, which looks like a dipole positioned 
at the two aforementioned in-phase sites. The difficulty 
in observing this evolution mode is explained by the fact 
that, in the course of the norm exchange, only ~ 0.1% 
of the total norm is actually transferred between the two 
sites. Furthermore, as mentioned earlier, the correspond- 
ing small real eigenvalue is completely suppressed by the 
surface (see panel (c) in Fig. [2]). It is worth noting that 
such stable three-site horseshoe surface structures may 
also be generated by the evolution of more complex un- 
stable waveforms, such as the five-site pyramids placed 
normally to the surface, see the bottom panel in Fig. [T2] 



V. CONCLUSIONS 

In this work, we have investigated localized modes in 
the vicinity of a two-dimensional surface, in the frame- 
work of the three-dimensional DNLS equation, which is 
a prototypical model of nonlinear dynamical lattices. We 
have found that the surface may readily stabilize local- 
ized structures that are unstable in the bulk (such as 
three-site horseshoes), and, on the other hand, it may in- 
hibit the formation of some other structures that exist in 
the bulk (such as vortices which are oriented normally to 
the surface, although ones parallel to the surface do exist 





t=182 * A ™ 



(c) t=0 




FIG. 12: (Color Online) The evolution of unstable pyramids. 
Panels (a), (b), and (c) display, respectively, the transforma- 
tion of a bulk pyramid, and of ones oriented normally and 
parallel to the surface, for e = 0.2. 



and have their stability region; a qualitative explanation 
to these features was proposed, based on the analysis 
of the interaction of the vortical state with its "mirror 
image"). The most typical surface- induced effect is the 
expansion of the stability intervals for various solutions 
that exist in the bulk and survive in the presence of the 
surface. This feature may be attributed to the decrease, 
near the surface, of the number of neighbors to which ex- 
cited sites couple, since the approach to the continuum 
limit, i.e., the strengthening of the linear couplings to the 
nearest neighbors, is responsible for the onset of the in- 
stability or disappearance of all the localized stationary 
states in the three-dimensional dynamical lattice. 

On the other hand, while the techniques elaborated in 
Rcfs. [HI, [13, f° r the analysis of localized states in 
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bulk lattices arc quite useful in understanding the domi- 
nant stability properties of the solutions, the surface gives 
rise to specific effects, such as the stabilization of higher- 
order solutions or the suppression of some types of vortex 



structures, which cannot be explained by these methods. 
Therefore, it would be very relevant to modify these tech- 
niques, which are based on the Lyapunov-Schmidt reduc- 
tions, so as to take the presence of the surface into regard. 
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